1 Введение в дисциплину. Моделирование при исследовании скважин и пластов
1.1 Тезисы для лекции
Связь исследований и моделирования через процесс принятия решений
Управление месторождениями строится на принятии решений. Принятие решений связано с необратимой тратой ресурсов с целью оптимизации некоторой целевой функции. Чаще всего целевая функция связана с экономикой проекта.
Исследования - инструмент для получения информации для принятия решений. Полезно разделить понятия данные и информация. Данные - то, что записано, информации - данные обработанные в контексте принятия решений.
Исследования основаны на моделях объекта управления. Модель - иструмент получения информации из данных. Также модель - инструмент принятия решений на основе информации. Модель - упрощенное представление реальности.
Модели используемые для принятия решений должны адекватно отражать реальность. Обеспечивается настройкой - адаптацией модели на известные данные и информацию. Такой процесс также можно считать исследованием. Почему?
Принятие решений на основе информации часто осуществляется на основе прогнозов - альтернативных вариантах принятия решений. Модель - инструмент для построения прогнозов.
Схема принятия решений с использованием моделей.
В узком смысле - исследование скважины, в частности гидродинамические исследования - соответствуют адаптации специализированной модели - решения уравнения фильтрации. Адаптация готовит модель для построения прогнозов и принятия решений.
Понимание всей цепочки создания ценности от данных до принятия решений оказывается полезным, например для обоснования необходимости проведения исследований и оценки ценности информации получаемой в ходе исследований
Кроме понимания всей цепочки важно развитие навыков проведения расчетов - адаптации различных моделей, построения прогнозов с использованием этих моделей. Этому будет посвящена значительная часть курса.
1.2 Инструменты для моделирования
Работа с моделями в ходе курса может проводиться с использованием различного ПО
python для построения аналитических моделей
Для работы можно использовать jupyter notebook или google colab или любой другой сервис для работы со скважинами. Потребуется освоить базовые библиотеки python - numpy, scipy, matplotlib, mpmath, pandas. Также будут использованы некоторые другие специализированные библиотеки.
Макросы MS Excel
Часть расчетов можно выполнить с использованием макросов MS Excel, например Unifloc VBA (https://github.com/unifloc/unifloc_vba)
Гидродинамический симулятор для построения численных моделей
Базовый вариант - тНавигатор. Версия начиная от 19. Можно будет использовать другие симуляторы - РН-КИМ, Schlumberger eclipse, OPM и другие.
Специализированное ПО для анализа ГДИС
Возможно, что то появится в ходе курса.
2 Уравнение фильтрации
Уравнение фильтрации многофазного потока в пористой среде - основная математическая модель пласта используемая в инженерных приложениях. Используется не само уравнение, а разнообразные частные решения, описывающие те или иные ситуации. Наиболее широко известным примером решения уравнения фильтрации, возможно является, численное решение которое строится трехмерными гидродинамическими симуляторами (eclipse, tnavigator и так далее).
Как правило, при выводе уравнения фильтрации используются следующие соотношения:
закон Дарси - определяет линейную зависимость скорости фильтрации от давления
\[
u_r=-\frac{k}{\mu}\frac{dp}{dr}
\tag{1}\]
уравнение неразрывности - или закон сохранения массы, показывает что в потоке не могут образовываться разрывы сплошности: \[
\frac{1}{r}\frac{\partial\left(r\rho u_r\right)}{\partial r}=-\varphi\frac{\partial p}{\partial t}
\tag{2}\]
уравнение состояния - задает зависимость плотности флюида от давления: \[
c_0=\frac{1}{\rho}\frac{\partial\rho}{\partial p}
\tag{3}\]
Опираясь на эти соотношения уравнение фильтрации в радиальной форме можно привести к виду для величин СИ.
\(c_t\) - общая сжимаемость породы и флюида, 1/атм
\(t\) - время, час
Вывод уравнения фильтрации основан на следующих предположениях:
пласт однородный и изотропный - пористость и проницаемость одинаковы во всем пласте и во всех направлениях и не зависят от давления
добывающая скважина вскрывает весь продуктивный горизонт и обеспечивает радиальный приток к скважине
пласт насыщен одним флюидом на всем протяжении
температура не меняется в пласте, изотермичность
2.1 Вывод уравнения фильтрации для линейного потока
Уравнение фильтрации основано на трех основных уравнениях:
законе сохранения массы
зависимости связывающей градиент давления с расходом, в простейшем виде это закон Дарси
уравнении состояния - зависимости связывающей сжимаемость флюида (и отчасти среды) и давлением
Выпишем эти уравнения, чтобы разобраться. Ориентируемся в обозначениях и подходе на [1]. Особое внимание на данном этапе стоит обратить на размерные коэффициенты в уравнениях. С одной стороны записывая уравнения в [СИ] их можно упростить и избежать части сложностей с переводными коэффициентами. С другой стороны во многих источниках переводные коэффициенты интесивно используются и это оказывается удобно при программной реализации алгоритмов. В данном описании ориентируясь на [1] приведем перевод коэффициентов в практические метрические единицы.
2.1.1 Закон сохранения массы
В дифференциальной форме, для потока в направлении \(x\)
\(f\) - размерный коэффициент, для практических метрических единиц \(f= 0.04167\), для американских промысловых \(f = 0.23394\)
\(\rho\) - плотность флюида, кг/м\(^3\), lbm/ft\(^3\)
\(q_x\) - объемный расход в направлении \(x\), м\(^3\)/сут, bbl/day, учитываем для практических метрических единиц \(1\) м\(^3\)/сут \(= 0.04167\) м\(^3\)/час и для промысловых американских \(1\) bbl/day \(= 0.23394\) ft\(^3\)/hr
x - координата, м, ft
\(\phi\) - пористость, доли единиц
\(t\) - время, часы
\(A\) - площадь, м\(^2\), ft\(^2\)
В перечне указаны размерности для промысловых метрических и американских промысловых единиц измерений, там где они отличаются.
Для практических метрических единиц (6) можно записать в виде \[
-0.04167 \dfrac{ \partial (\rho q_x) }{\partial x} = A \dfrac{\partial(\rho \phi)}{\partial t}
\tag{7}\]
2.1.2 Закон Дарси
Простейшим вариантом задания зависимости между измением давления и потоком флюида является линейная зависимость - закон Дарси. При этом линейный множитель в зависимости - функция как коллектора (проницаемость), так и свойств флюида (вязкость).
Закон Дарси в направлении \(x\). \[
q_x = - \dfrac{k_x A}{f\mu} \dfrac{\partial p}{\partial x}
\tag{8}\]
где
\(f\) - размерный коэффициент, для практических метрических единиц \(f= 115.74\), для американских промысловых \(f = 887.2\)
\(k_x\) - проницаемость в направлении движения потока, мД
\(A\) - площадь, м\(^2\), ft\(^2\)
\(\mu\) - динамическая вязкость, сП
\(p\) - давление, атм, psi
\(x\) - координата движения потока, м
В практических метрических единицах измерения Уравнение 8 будет иметь вид \[
q_x = - \dfrac{k_x A}{115.74\mu} \dfrac{\partial p}{\partial x}
\tag{9}\]
Пытаясь проверить корректность вывода переводных коэффициентов необходимо учесть что для единиц измерение СИ получим \(f=1\), а необходимые переводные коэффициенты можно записать как:
Принцип неразрывности показывает, что для определенного объема пласта, масса флюида которое втекла в контрольный объем пласта минус масса которая вытекла равна массе которая накопилась в объеме.
Сохранение импульса или закон Дарси можно выразить соотношением
\[
u_r=-\frac{k}{\mu}\frac{d^p}{dr}
\tag{15}\]
Закон Дарси здесь используется как псевдоустановившаяся аппроксимация обобщенного уравнения сохранения импульса (то есть слагаемым отвечающим за накопление импульса пренебрегаем). Это предположение справедливо если не учитывать возмущения давления в среде двигающиеся со скоростью звука. Все изменения давления описываемые моделью связаны с локальными изменениями градиента давления за счет ламинарного режима потока. Хотя далее в модели будет учтена сжимаемость системы всеми “звуковыми” эффектами в системе мы пренебрегаем.
Поток предполагается горизонтальным, поэтому давление \(p\) может быть использовано в качестве потенциала потока, гравитационными силами пренебрегаем.
Уравнение 16 – дифференциальное уравнение в частных производных описывающее нестационарный поток однофазного флюида в пористой среде при ламинарном потоке.
Вообще говоря приведенное уравнение является нелинейным, так как плотность \(\rho = \rho(p)\) и вязкость \(\mu = \mu(p)\) являются функциями давления. Уравнение содержит две зависимые переменные - давление \(p\) и плотность \(\rho\). Поэтому для его решения необходимо задать еще одно соотношение, каковым может быть уравнение состояния флюида связывающее плотность флюида и давление \(\rho = \rho(p)\).
2.2.2 Флюид постоянной сжимаемости
Нестационарное поведение давления в пласте связано со сжимаемостью системы. При изменении давления в какой то точке, часть флюида сжимается, происходит накопление или отдача флюида, что вызывает задержку в распространении изменения флюида. Несмотря на то, что сжимаемость флюидов и породы малы и во многих случаях ими можно пренебречь, это не верно для пластовых систем для добычи нефти. Большие объемы пласта и флюидов и высокие давления компенсируют малость сжимаемости и требуют ее учета.
Для однофазного флюида разумным является предположение постоянства сжимаемости.
хотя пористость здесь является функцией давления - в первом приближении мы можем считать ее константой равной пористости при некотором среднем давлении в пласте. Это справедливо для маленькой сжимаемости породы, что верно почти всегда.
Уравнение для сжимаемости можно еще уточнить, учтя что в пласте могут находится различные флюиды - вода и нефть с насыщенностями \(s_w\) и \(s_o\)
тогда
\[
c_l = s_o c_o + s_w c_w
\]
тогда можно ввести общую сжимаемость системы
\[
c_t = c_l + c_f = s_o c_o + s_w c_w + c_f
\]
Заметим, что проницаемость k в законе Дарси это не абсолютная проницаемость, но относительная проницаемость по нефти при насыщенности водой соответствующей связанной воде.
\[
k= k_o (s_{wc})
\]
2.2.4 Линеаризация уравнения фильтрации
Раскрыв производную в левой части Уравнение 17 и предположив, что \(\dfrac{\partial p}{\partial r}\) мало а следовательно слагаемым \(r \rho c_t \left( \dfrac{\partial p}{\partial r} \right)^2\) можно пренебречь, получим линеаризованное уравнения фильтрации
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\)
Для практических расчетов удобнее бывает использовать безразмерные переменные полученные для практических метрических единиц измерения. \[
r_D = \frac{r}{r_w}
\]
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния \(p_D(r_D, t_D)\)
2.3.1 Расчет безразмерных переменных в Unifloc VBA
Несмотря на простоту определений безразмерных переменных их часто приходится применять при проведении расчетов. Поэтому в надстройке Unifloc VBA реализован набор функций расчета безразмерных переменных.
Описания функций и из аргументов можно найти в руководстве пользователя Unifloc VBA
3 Стационарные решения уравнения фильтрации
Широкое распространение на практике получили стационарные решения уравнения фильтрации. Приведем некоторые из них.
3.1 Решение для постоянного давления на круговой границе
Рассматривается самая простая модель работы добывающей скважины - радиальная стационарная фильтрация в однородном изотропном пласте круговой формы. Скважина находится в центре пласта (Рисунок 1). На границе пласта поддерживается постоянное давление. Фактически это означает, что через границу пласта идет поток жидкости, уравновешивающий дебит скважины.
Решение можно получить как решения уравнения фильтрации, учитывая стационарность потока
Рисунок 1: Схема радиального притока к скважине при наличии постоянного давления на границе
Или можно получить его непосредственно из закона Дарси, который должен быть приведен к радиальной форме и в таком варианте известен как Формула Дюпюи. Для приведенной конфигурации можно записать закон Дарси в форме.
\(q\) - объемные дебит скважины в рабочих условиях, м\(^3\)/сут
\(r\) - радиус - расстояние от центра скважины, м
\(r_e\) - радиус зоны дренирования, на котором поддерживается постоянное давление, м
\(r_w\) - радиус скважины, на котором замеряется забойное давление, м
\(P\) - давление, атм
\(P_e\) - давление на внешнем контуре дренирования, атм
\(P_w\) - давление на забое скважины, атм
\(k\) - проницаемость, мД
\(\mu\) - вязкость нефти в зоне дренирования, сП
Далее если не указано особо будем использовать практические метрические единицы.
3.1.1 Учет скин-фактора
Скин-фактор — гидродинамический параметр, характеризующий дополнительное фильтрационное сопротивление течению флюидов в околоскважинной зоне пласта, приводящее к изменению добычи (дебита) по сравнению с совершенной (идеальной) скважиной. Скин-фактор может приводить как к снижению дебита (например при загрязнении ПЗС), так и увеличению (образование высокопроводящих каналов в ПЗС).
Концепция скин-фактора получила широкое распространение на практике. Все инженеры-нефтяники знают этот параметр и оперируют им на практике.
Изначально скин-фактор был введен как параметр учитывающий изменение проницаемости (загрязнение) призабойной зоны при расчете производительности скважины. Такое загрязнение может быть вызвано различными причинами:
проникновением бурового раствора в пласт и блокировкой поровых каналов;
набуханием глин при контакте с фильтратом бурового раствора;
химическим осаждением элементов бурового раствора, жидкости глушения или пластовых флюидов в призабойной зоне скважины, например осаждением солей или асфальтенов;
продвижением песчаных частиц к стволу скважины;
повреждением породы при перфорации;
другими причинами.
Для модели загрязненной призабойной зоны величину скин-фактора можно выразить формулой Хокинса . Скин-фактор для плоскорадиального установившегося потока несжимаемой жидкости:
\[
S =\left( \frac{k}{k_s} -1\right)\ ln\frac{r_s}{r_w}
\tag{22}\]
здесь:
\(k_s\) - проницаемость в загрязненной ПЗП;
\(k\) - однородная проницаемость по всему пласту;
\(r_s\) - радиус загрязненной зоны;
\(r_w\) - радиус скважины.
Концепция скин-фактора оказалась удобной для описания характеристики соединения скважины и пласта и была распространена на другие случаи, когда производительность скважины могла отличаться от производительности идеальной скважины:
для горизонтальных скважин;
для скважин вскрывающих пласт под углом;
для скважин пересеченных трещиной ГРП;
для скважин вскрытых перфорацией и учета гидравлического сопротивления потока на перфорационных отверстиях;
другими причинами.
Для многих подобных случаев предположение о радиальном притоке к скважине не верно, но величину скин-фактора используют, так как она позволяет сравнить производительность скважины со сложным заканчиванием с простой вертикальной скважиной. В таких случая говорят о псевдорадиальном скин-факторе - такой величине скин-фактора \(S\), которая обеспечила бы такую же производительность для вертикальной скважины полностью вскрывающей пласт.
Для стационарной радиальной модели притока учет скин-фактора приведен к следующим соотношениям:
3.1.2 Решение для постоянного давления на круговой границе с учетом среднего давления в области дренирования
Рисунок 2: Среднее давление в области дренирования
В приведенное выражение входит значение давления на контуре, которым не всегда бывает удобно пользоваться. В практических случаях значение на контуре трудно оценить, контур зоны дренирования может значительно отличаться от кругового, да и радиус оценить может быть сложно. Удобнее пользоваться средним давлением в зоне дренирования \(\bar{P}\), которое может быть оценено по материальному балансу (смотри рисунок \(\ref{ris:radial_inflow_steady_state_average_pressure}\)).
Вывод уравнения фильтрации для постоянного давления на границе с использованием среднего давления
В этом случае выражение для дебита примет вид \[
q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{1}{2}+ S \right)}
\]
3.2 Решение для круговой непроницаемой границы
Схема модели радиального притока для условия непротекания на круговой границе приведена на Рисунок 3
Рисунок 3: Схема радиального притока к скважине при наличии непроницаемой границы
При условии непротекания давления на границе условия стационарности (неизменности давления) не достигаются. При работе скважины с постоянным дебитом забойное давление будет постоянно снижаться. Однако начиная с некоторого момента, когда влияние скважины достигнет границ - давление в всей области дренирования начнет снижаться равномерно (смотри Рисунок 4).
Рисунок 4: Изменение давления на границе и на забое скважины во времени
Такой режим, при котором забойное давление меняется, но перепад давления \(P_e - P_w\) остается постоянным называют псевдо-установившимся режимом работы (pss - pseudo steady state).
Для псевдо-установившегося режима можно записать выражение
\[
q=\frac{kh\left(P_e-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{1}{2} + S\right)}
\] где
\(q\) - объемные дебит скважины в рабочих условиях, м3/сут;
\(r\) - радиус - расстояние от центра скважины, м;
\(r_e\) - радиус зоны дренирования, на котором поддерживается постоянное давление, м;
\(r_w\) - радиус скважины, на котором замеряется забойное давление, м;
\(P\) - давление, атм;
\(P_e\) - давление на внешнем контуре дренирования, атм;
\(P_w\) - давление на забое скважины, атм;
\(k\) - проницаемость, мД;
\(\mu\) - вязкость нефти в зоне дренирования, сП.
3.2.1 Решение для круговой непроницаемой границы с учетом среднего давления в зоне дренирования
Аналогично случаю для постоянного давления на границе можно переписать выражение с использованием среднего давления в области дренирования.
\[
q=\frac{kh\left( \bar{P}-P_w\right)}{ 18.41 \mu\left(\ln{\dfrac{r_e}{r_w}} - \dfrac{3}{4}+ S \right)}
\]
Вывод уравнений для псевдо-установившегося режима работы
3.2.2 Стационарные решения для вертикальной скважины в резервуаре произвольной формы
Здесь уравнения и методы расчета для горизонтальных, наклонно направленных скважин, скважин с ГРП, горизонтальных скважин с МГРП.
Простое решение для задачи стационарного притока к вертикальной скважине в однородном изотропном пласте круговой формы с постоянным давлением на границе имеет вид
который удобен для расчёта распределения давления в пласте \(P_r\) на произвольном расстоянии от скважины \(r\). В выражении (2) задано граничное значение давления \(p_e\) на контуре \(r_e\). Расчёт позволит найти любое значение внутри контура, в том числе и забойное давление \(P_{wf}\) на \(r=r_w\)
Выражение можно переписать \[
P_{r} = P_{wf} + 18.41\dfrac{ Q\mu B }{kh} \left[ \ln\dfrac{r}{r_w} +S \right]
\]
где по известному дебиту и забойному давлению можно найти давление в пласте. При известном пластовом давлении можно оценить радиус контура на котором оно достигается.
3.4 Расчет на python
Расчет на python может быть реализован на основе библиотек numpy scipy
Show the code
"""Импортируем библиотеки для расчетов. numpy используем для работы с массивами и подготовки данных для построения графиков.matplotlib используем для построения графиковscipy для решения линейных уравнений"""import numpy as npimport matplotlib.pyplot as pltimport scipy
Для удобства дальнейшего изложения и использования расчетных функций при создании функций и переменных на языке python названия формируются по следующим принципам:
сначала указывается, что расчитывается в функции, в данном случае - давление \(p\)
потом указываются пояснения - в данном случае p_ss - steady state pressure
в конце указывается размерность в которой ожидается получаение ответа - в данном случае atma - абсолютные атмосферы.
Show the code
"""Определим функции для расчета стационарного решения"""def dp_ss_atm(q_liq_sm3day =50, mu_cP =1, b_m3m3 =1.2, kh_mDm =40, r_e_m =240, r_m =0.1):""" функция расчета перепада давления в произвольной точке пласта на расстоянии r_m от центра скважины для стационарного решения - q_liq_sm3day - дебит жидкости на поверхности в стандартных условиях - mu_cP - вязкость нефти (в пластовых условиях) - B_m3m3 - объемный коэффициент нефти - kh_mDm - kh пласта - r_e_m - радиус контрура питания, м - r_m - расстояние на котором проводится расчет, м """return18.42* q_liq_sm3day * mu_cP * b_m3m3/ kh_mDm * np.log(r_e_m/r_m)def p_ss_atma(p_res_atma =250, q_liq_sm3day =50, mu_cP =1, b_m3m3 =1.2, k_mD =40, h_m =10, r_e_m =240, r_m =0.1):""" функция расчета давления в произвольной точке пласта на расстоянии r_m от центра скважины для стационарного решения - p_res_atma - пластовое давление, давление на контуре питания - q_liq_sm3day - дебит жидкости на поверхности в стандартных условиях - mu_cP - вязкость нефти (в пластовых условиях) - B_m3m3 - объемный коэффициент нефти - k_mD - проницаемость пласта - h_m - мощность пласта, м - r_e_m - радиус контрура питания, м - r_m - расстояние на котором проводится расчет, м """return p_res_atma - dp_ss_atm(q_liq_sm3day = q_liq_sm3day, mu_cP = mu_cP, b_m3m3 = b_m3m3, kh_mDm = k_mD * h_m, r_e_m = r_e_m, r_m = r_m)
Функции расчетов могут быть использованы для построения графиков, например с использованием matplotlib
Show the code
"""Построим график распределения давления в пласте"""# формируем массив расстояний для которых будем проводить расчетr_arr = np.logspace(np.log10(0.1), np.log10(240), 500) # рассчитываем массив давлений на соответствующих расстояниях# для расчета используется векторный расчет numpy # для примера показана передача всех аргументов созданной функцииp_arr = p_ss_atma(p_res_atma =250, q_liq_sm3day =50, mu_cP =1, b_m3m3 =1.2, k_mD =40, h_m =10, r_e_m =240, r_m = r_arr)fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))# рисуем график в обычных координатахax1.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивамax1.plot(-r_arr, p_arr) # настраиваем графикax1.set_title('Воронка депрессии - распределение давления в пласте')ax1.set_xlabel('r, m')ax1.set_ylabel('p, atma')# рисуем график в логарифмических координатахax2.plot(r_arr, p_arr) # команда отрисовки графика по заданным массивамax2.plot(-r_arr, p_arr) # настраиваем графикax2.set_title('Воронка депрессии - распределение давления в пласте')ax2.set_xlabel('r, m')ax2.set_ylabel('p, atma')ax2.set_xscale('symlog', linthresh=0.1, linscale=0.5)plt.show()
Рисунок 6: Распределение давления в круговом пласте
3.4.1 Формула Дюпюи в декартовых координатах
Show the code
"""Построим карту или сетку распределения давления в пластес использованием векторных функций numpy"""# зададим параметры воронки депрессииpres =250r_e =300# зададим координатную сетку основываясь на параметрахx = np.linspace(-r_e, r_e, 300)y = np.linspace(-r_e, r_e, 300)# рассчитаем вспомогательные вектора для построения сеткиxv, yv = np.meshgrid(x, y)# зададим координаты скважиныxwell1 =0ywell1 =0# рассчитаем значение давлений во всех точках сетки# расчет ведется для матрицы координата с использованием векторных возможностей numpyp_mesh = p_ss_atma(r_m=((xv-xwell1)**2+ (yv-ywell1)**2)**0.5, p_res_atma=pres, r_e_m=r_e)# удалим значения за контуром, так как в данном случае они не имеют смыслаp_mesh[np.where(p_mesh > pres)] = pres# построим отображение в виде контурной картыplt.rcParams["figure.figsize"] = (7,7)plt.contour(x, y, p_mesh, levels=100)plt.show()
Рисунок 7: Карта изолиний давления в пласте
Для построения карты распределения давлений в пласте полезно вспомнить, что расстояние от скважины с координатами \((x_{well}, y_{well})\) до произвольной точки пласта с координатами \((x,y)\) можно найти по формуле \[
r=\sqrt{ (x-x_{well})^2 + (y-y_{well})^2 }
\]
Тогда выражение для расчета давления в любой точке пласта примет вид
Простой вариант расчета - можно создать пустую матрицу со значениями давления по сетке и перебирая все точки на сетке/матрице рассчитать давления
3.5 Суперпозиция для стационарного решения
3.5.1 Суперпозиция для нескольких скважин с постоянным дебитом
Show the code
"""Построим карту или сетку распределения давления в пластедля ряда/галереи скважин"""# зададим параметры воронки депрессииpres =250r_e =300# зададим координаты скважины (всего будет 5 скважин у которых меняется только х координата)xwell = [-100, -50, 0 , 50, 100]ywell =0# зададим координатную сетку основываясь на параметрахx = np.linspace(-r_e, r_e, 300)y = np.linspace(-r_e, r_e, 300)# рассчитаем вспомогательные вектора для построения сеткиxv, yv = np.meshgrid(x, y)# зададим пустой список матриц с перепадами давлений от каждой отдельной скважиныp_mesh_i=[]# для каждой скважины найдем ее влияние на давленияfor xi in xwell:# рассчитаем матрицу расстояний от элементов сетки до скважины i r_well = ((xv-xi)**2+ (yv-ywell)**2)**0.5# рассчитаем значение давлений во всех точках сетки для скважины i p_mesh_i_ = p_ss_atma(r_m=r_well, p_res_atma=pres, r_e_m=r_e)# удалим значения за контуром, так как в данном случае они не имеют смысла p_mesh_i_[np.where(p_mesh_i_ > pres)] = pres# сохраним влияние скважины i в списке матриц влияния отдельных скважин p_mesh_i.append(p_mesh_i_)# найдем суммарное влияние все скважинp_mesh_sum =0for p_mesh_i_ in p_mesh_i:# найдем сумму решений, помним что суммировать можно депрессии p_mesh_sum = pres - p_mesh_i_ + p_mesh_sump_mesh_sum = pres - p_mesh_sumfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))# построим отображение в виде контурной картыax1.contourf(x, y, p_mesh_sum, levels=100)# построим отображение в виде контурной картыax2.contour(x, y, p_mesh_sum, levels=15)plt.show()
Рисунок 8: Карта изолиний давления в пласте для нескольких скважин с постоянным дебитом
Show the code
"""Построим график изменения давления по линии расположения скважин"""xline = np.linspace(-300, 300, 10000)yline =0pline = np.zeros_like(xline)for xi in xwell: r_well = ((xline-xi)**2+ (yline-ywell)**2)**0.5 r_well = np.where(r_well <0.1, 0.1, r_well) pline = pline + p_ss_atma(r_m=r_well, p_res_atma=0, r_e_m=r_e, q_liq_sm3day =50, mu_cP =1, b_m3m3 =1.2, k_mD =40, h_m =10)pline = pres + plineplt.plot(xline, pline)plt.show()
Рисунок 9: Распределение давления в пласте для нескольких скважин с постоянным дебитом по линии скважин
Для стационарного решения работает принцип суперпозиции - сумма двух решений также будет решением, это позволяет построить карту для нескольких скважин. Давление в любой точке пласта можно найти по формуле
Если считать забойные давления \(P_{wf.j}\) известными а дебиты скважин \(Q_i\) не известными, тогда выражение (6) можно рассматривать как систему линейных алгебраических уравнений вида
\[
AX = B
\]
Где \[
A_{[i,j]} = 18.41\dfrac{ \mu B }{kh} \left[ \ln\dfrac{r_e}{\sqrt{ (x_{w.j}-x_{w.i})^2 + (y_{w.j}-y_{w.i})^2 }} +S \right]
\]
\[
B_{[j]}=P_{res} - P_{wf.j}
\]
такую систему можно решить например с использованием пакета scipy.linalg
Show the code
# зададим параметры воронки депрессииpres =250r_e =300# зададим координатную сетку основываясь на параметрахx = np.linspace(-r_e, r_e, 500)y = np.linspace(-r_e, r_e, 500)# рассчитаем вспомогательные вектора для построения сеткиxv, yv = np.meshgrid(x, y)# зададим координаты скважины (всего будет 5 скважин у которых меняется только х координата)xwell = [-100, -50, 0 , 50, 100]ywell = [0,0,0,0,0]pwell = [100, 100, 100, 100, 100] # забойные давления скважинqwell = [0,0,0,0,0] # дебиты скважин, должны быть найдены# создадим заготовки матрицы А и вектора В по формулам (4) и (5)A = np.zeros((5,5))B = np.zeros(5)# сформируем матрицу Аfor i,v inenumerate(xwell):for j,v inenumerate(ywell):# найдем расстояния от одной скважины до другой# ищем расстояния между центрами скважин r_ij = ((xwell[i]-xwell[j])**2+ (ywell[i]-ywell[j])**2)**0.5# если расстояние 0 значит ищем влияние скважины саму на себя# тогда подставляем радиус скважиныif r_ij ==0: r_ij =0.1# чтобы воспользоваться ранее заданной формулой в виде функции# вызовем ее с нулем пластовым давлением и единичным дебитом A[i,j] =- p_ss_atma(p_res_atma =0, q_liq_sm3day =1, mu_cP =1, b_m3m3 =1.2, k_mD =40, h_m =10, r_e_m =240, r_m = r_ij) B[i] = pres - pwell[i]# найдем решение = значения дебитов при которых забойные равны заданным значениямqwell = scipy.linalg.solve(A,B)# напечатаем найденные дебитыprint(qwell)# зададим пустой список матриц с перепадами давлений от каждой отдельной скважиныp_mesh_i=[]# для каждой скважины найдем ее влияние на давленияfor i,xi inenumerate(xwell):# рассчитаем матрицу расстояний от элементов сетки до скважины i r_well = ((xv-xwell[i])**2+ (yv-ywell[i])**2)**0.5# для красоты отрисовки карты давлений пренебрежем значениями в радиусе одного метра от скважин r_well[r_well<1]=0.1# рассчитаем значение давлений во всех точках сетки для скважины i p_mesh_i_ = p_ss_atma(r_m=r_well, p_res_atma=pres, r_e_m=r_e, q_liq_sm3day=qwell[i])# удалим значения за контуром, так как в данном случае они не имеют смысла p_mesh_i_[np.where(p_mesh_i_ > pres)] = pres# сохраним влияние скважины i в списке матриц влияния отдельных скважин p_mesh_i.append(p_mesh_i_)# найдем суммарное влияние все скважинp_mesh_sum_new =0for p_mesh_i_ in p_mesh_i:# найдем сумму решений, помним что суммировать можно депрессии p_mesh_sum_new = pres - p_mesh_i_ + p_mesh_sum_newp_mesh_sum_new = pres - p_mesh_sum_newfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))# построим отображение в виде контурной картыax1.contourf(x, y, p_mesh_sum_new, levels=100)# построим отображение в виде контурной картыax2.contour(x, y, p_mesh_sum_new, levels=15)plt.show()
Рисунок 11: Распределение давления в пласте для нескольких скважин с постоянным забойным давлением по линии скважин
3.5.3 Расчета поля давлений для нескольких скважин на стационарном режиме с использованием классов
Классы в python помогают задавать сложные структуры данных (инкапсуляция) и работать с ними как с едиными объектами. Скважины с их параметрами, а также группы скважин могут быть заданы с использованием классов.
Класс для расчета влияния одной скважины - объединяет все необходимые данные для скважины - позволяет рассчитать депрессию и давления в любой точке пласта с учетом радиуса влияния и радиуса скважины
Show the code
class Well:""" класс Well описывает влияние одной скважины на поле давления используется стационарное решение """def__init__(self, name, xw_m, yw_m) ->None:self.x_m = xw_mself.y_m = yw_mself.name = nameself.kh_mDm =10.0self.b_m3m3 =1.0self.mu_cP =1.0self.skin =0self.r_e_m =250.0self.r_w_m =0.1self.q_liq_sm3day =10def calc_dp_ss_atm(self, x_m, y_m):""" расчет перепада давления в произвольной точке пласта с заданными координатами для стационарного решения уравнения фильтрации x_m, y_m координаты в которых рассчитывается давление """ r_m = np.sqrt( (self.x_m - x_m)**2+ (self.y_m - y_m)**2 ) r_m = np.where(r_m >=self.r_e_m, self.r_e_m ,r_m) r_m = np.where(r_m <self.r_w_m, self.r_w_m ,r_m)return dp_ss_atm(q_liq_sm3day=self.q_liq_sm3day, mu_cP=self.mu_cP, b_m3m3=self.b_m3m3, kh_mDm=self.kh_mDm, r_e_m=self.r_e_m, r_m=r_m)def calc_dp_ss_well_atm(self, w):returnself.calc_dp_ss_atm(w.x_m, w.y_m)def calc_p_ss_atma(self, x_m, y_m, p_res_atma):""" функция расчета давления в произвольной точке пласта с заданными координатами для стационарного решения уравнения фильтрации x_m, y_m координаты в которых рассчитывается давление """return p_res_atma -self.calc_dp_ss_atm(x_m, y_m)class Wells:""" класс для группы взамодействующих скважин """def__init__(self, name:list, xw_m:list, yw_m:list, kh_mDm=10, b_m3m3=1, mu_cP=1, r_e_m =250, r_w_m=0.1) ->None:self.well_dict: dict[str, Well] = {} for (n, x, y) inzip(name, xw_m, yw_m):self.well_dict.update({n: Well(n,x,y)})self.set_kh_mDm(kh_mDm)self.set_pvt(b_m3m3=b_m3m3, mu_cP=mu_cP)self.set_r_e_m(r_e_m)self.set_r_w_m(r_w_m)self.p_res_atma =250def set_qliq_sm3day(self, q_liq_list_sm3day:list):for (wname, q) inzip(self.well_dict, q_liq_list_sm3day):self.well_dict[wname].q_liq_sm3day = qdef set_kh_mDm(self, kh_mDm):for wn inself.well_dict:self.well_dict[wn].kh_mDm = kh_mDmdef set_r_e_m(self, r_e_m):for wn inself.well_dict:self.well_dict[wn].r_e_m = r_e_mdef set_r_w_m(self, r_w_m):for wn inself.well_dict:self.well_dict[wn].r_w_m = r_w_mdef set_pvt(self, mu_cP, b_m3m3):for wn inself.well_dict:self.well_dict[wn].b_m3m3 = b_m3m3self.well_dict[wn].u_cP = mu_cPdef calc_p_ss_atma(self, x_m, y_m): dp =0for wn inself.well_dict: dp = dp +self.well_dict[wn].calc_dp_ss_atm(x_m, y_m)returnself.p_res_atma - dpdef estimate_qliq_from_pwf(self, wn_list:list, pwf_list:list): dim =len(wn_list) # размерность матрицы для расчета дебитов# создадим заготовки матрицы А и вектора В по формулам (4) и (5) A = np.zeros((dim,dim)) B = np.zeros(dim)for i, wni inenumerate(wn_list):self.well_dict[wni].q_liq_sm3day =1for j, wnj inenumerate(wn_list): A[i,j] =self.well_dict[wni].calc_dp_ss_well_atm(self.well_dict[wnj]) B[i] =self.p_res_atma - pwf_list[i] q_list = scipy.linalg.solve(A,B)for i, wni inenumerate(wn_list):self.well_dict[wni].q_liq_sm3day = q_list[i]xlist = [300, -300, 50]ylist = [0, 0, 200]names = ['1', '2', '3']qlist = [1, -20, 10]ww = Wells(names, xlist, ylist, r_e_m=500)ww.set_qliq_sm3day(qlist)xp = np.linspace(-500, 500, 500)yp = np.linspace(-500, 500, 500)xv, yv = np.meshgrid(xp, yp)p_mesh = ww.calc_p_ss_atma(xv, yv)fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))# построим отображение в виде контурной картыax1.contour(xp, yp, p_mesh, levels=50)for x,y,n inzip(xlist, ylist, names): ax1.text(x+10, y+10, n)ax1.scatter(xlist, ylist, s =100)ax2.streamplot(xp, yp, -np.gradient(p_mesh,axis=1), -np.gradient(p_mesh,axis=0), density=2)for x,y,n inzip(xlist, ylist, names): ax2.text(x+10, y+10, n)ax2.scatter(xlist, ylist, s =100)plt.show()
Рисунок 12: Расчета поля давлений для нескольких скважин с постоянным дебитом
3.5.4 Расчета поля давлений для нескольких скважин на стационарном режиме для заданных забойных давлений
Уравнение 25 можно интерпретировать следующим образом. Параметр \(T\) отвечает за свойства пласта и флюида на которые трудно повлиять в ходе эксплуатации. Это то, что дала природа в точке где находится скважина. Депрессия \(\Delta P\) – параметр которым можно управлять в ходе эксплуатации регулируя забойное давление. Например за счет установки насоса и задания параметров его работы. На этом параметре должно быть сосредоточено основное внимание при анализе работы скважины. Параметр \(J_D\) – определяет качество соединения скважины с пластом или качество заканчивания. Его мы можем выбирать при строительстве скважины и можем менять в ходе эксплуатации проводя ГТМ, хотя и достаточно большой ценой. Поскольку мы можем влиять на \(J_D\) важно понимать, какое оптимальное значение продуктивности можно достичь на конкретной скважине и как его можно изменить.
Задачей гидродинамических исследований является установление величин \(T\) и \(J_D\), хотя традиционно речь ведется об определении проницаемости \(k\) и скин-фактора \(S\).
3.7 Задания для самостоятельной работы
Для совершенствования навыков работы с python выполните следующие задания:
Постройте график распределения давления в пласте для композитного пласта. В композитном пласте на расстоянии \(r<r_1\) проницаемость равна \(k=k_1\), а для \(r>=r_1\), \(k=k_2\).
Постройте двумерную тепловую или контурную карту распределения давления в пласте для моделей однородного и композитного пласта.
Рассчитайте среднюю величину давления в круговой области дренирования для однородного пласта. Насколько среднее давление в круговой области дренирования будет отличаться от давления на контуре. Чему будет равен коэффициент \(S\) в выражении \(Q=\dfrac{kh}{18.41\mu B} \dfrac{P_{res}-P_{wf}}{ln(\dfrac{r_e}{r_w})+S}\) при использовании вместо давления на контуре среднего давления? Постройте график, на котором будет отображаться распределение давления в зоне дренирования и величина среднего давления (в виде линии).
Для примера с несколькими скважинами имитирующими трещину ГРП рассчитайте дебиты скважин таким образом, чтобы забойное давление на всех скважинах было одинаковым. Постройте графики распределения давления в пласте. Постройте график дебитов вдоль “скважины”.
4 Нестационарные решения
Для установившегося режима фильтрации давление в пласте не меняется. Для псевдо-установившегося режима постоянным остается перепад давления между пластом и забоем. После запуска, остановки или изменения режима работы скважины эти условия не выполняются. Давление в различных точках пласта может меняться по разному. Такой режим называют неустановившимся, а решения его описывающие нестационарными (зависят от времени).
Неустановившиеся решения уравнения фильтрации (transient solutions) представляют значительный практический интерес во многих задачах, включая задачи интерпретации ГДИС. В тоже время они относительно сложны и требуют применения компьютерных алгоритмов. В данном пособие проведение расчетов иллюстрируется с использованием python и макросов для Excel – Unifloc VBA.
4.1 Решение линейного стока
Для решения уравнения фильтрации - линейного дифференциального уравнения в частных производных второго порядка необходимо задать начальные и граничные условия.
Самое простое решение можно получить для случая вертикальной скважины бесконечно малого радиуса запускающейся с постоянным дебитом. Условия соответствующие этому случаю можно выразить следующим образом:
Начальное условие. До запуска скважины в момент времени \(t_D = 0\) давление в пласте равно начальному во всех точках \(p=p_i\)\[
t_D < 0, p_D = 0
\]
Граничное условие на скважине. Условие постоянства дебита на скважине можно трансформировать в граничное условие опираясь на закон Дарси. \[
\lim_{r_D \to 0} {r_D \frac{\partial p_D}{\partial r_D}} = -1
\]
Граничное условие на бесконечном расстоянии от скважины. Давление в пласте на бесконечно большом расстоянии от скважины равно начальному. \[
r_D = \infty, p_D = 0
\]
В этом случае решение может быть выражено через функцию интегральной экспоненты \[
p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right)
\]
где \(-Ei(-x)\) - интегральная показательная функция, рисунок \(\ref{ris:Ei_plot_1}\).
# импортируем библиотеки для расчетов# numpy используем для работы с массивами и подготовки данных для построения графиков. # Также в некоторых функциях используем возможности векторных расчетов numpyimport numpy as np# matplotlib используем для построения графиковimport matplotlib.pyplot as plt# scipy.special используем как альтернативный вариант расчета специальных функцийimport scipy.special as sc# для скорости и удобства используем sc.expi# построим ветки функции для положительных и отрицательных аргументов раздельноx = np.arange(1e-5,3,0.01)x1 = np.arange(-3,-1e-5,0.01)plt.plot(x, sc.expi(x))plt.plot(x1, sc.expi(x1))plt.title("График функции Ei")plt.xlabel("x")plt.ylabel("Ei(x)")plt.show()
Рисунок 14: График функции интегральной экспоненты Ei
Часто для проведения расчетов, особенно с использованием компьютерных библиотеке расчетов, бывает удобнее пользоваться модифицированной интегральной показательной функцией \(Ei_1(x)\) или \(E_1(x)\) или \(Ei_n(x)\) при \(n=1\). \[
Ei_n(x) = \int\limits_{1}^{\infty}\frac{e^{-tx}}{t^n}\,\mathrm dt
\]
График интегральной показательной функции \(Ei_1(x)\) приведен на Рисунок 14. Для вещественных положительных \(x\in\mathbb R, x>0\) верно \(E_1(x) = - Ei( -x)\)
Функцию интегральной экспоненты можно представить в виде ряда.
# зададим логарифмическое распределение точек вблизи нуля для построения графикаx = np.logspace(-10, 1, 100)plt.rcParams["figure.figsize"] = (8,3)fig, (ax1, ax2) = plt.subplots(1, 2)ax1.plot(x, sc.expi(-x))ax1.plot(x, np.log(x) +0.57721566481)#ax1.set_title('Сравнение Ei и ln для обычных координат')ax2.plot(x, sc.expi(-x))ax2.plot(x, np.log(x) +0.57721566481)ax2.set_xscale('log')#ax2.set_title('Сравнение Ei и ln для полулогарифмических координат')plt.show()
Рисунок 15: Сравнение Ei и ln для обычных и полулогарифмических координат
Из приведенного выражения можно сделать выводы, что для маленьких значений аргумента функция интегральной экспоненты \(E_1(x)\) может быть аппроксимирована логарифмической зависимостью.
\[
E_1(x) = -ln(x) - \gamma
\]
График сравнения функций \(E_1(x)\) и \(ln(x)\) показан на Рисунок 15. Видно, что хорошей аппроксимация будет только для маленьких значений аргумента \(x < 0.01\). Но для решения уравнения фильтрации именно эта зона представляет наибольший интерес.
Представление интегральной экспоненты в виде логарифмической аппроксимации удобно на практике, так как логарифм легче вычислять. В большинстве языков программирования и инструментов для проведения расчетов расчет логарифма реализован по умолчанию. А для расчета интегральной экспоненты, часто приходится предпринимать дополнительные шаги.
Решение уравнения фильтрации для линейного стока с учетом логарифмической аппроксимации можно представить в виде
Решение этого уравнения - функция безразмерного давления от безразмерных времени и расстояния $p_D(r_D, t_D) $
Для решения уравнения фильтрации - линейного дифференциального уравнения в частных производных второго порядка необходимо задать начальные и граничные условия. Самое простое решение можно получить для случая вертикальной скважины бесконечно малого радиуса запускающейся с постоянным дебитом. Условия соответствующие этому случаю можно выразить следующим образом
начальное условие. До запуска скважины в момент времени \(t_D = 0\) давление в пласте равно начальному во всех точках \(p=p_i\)\[
t_D < 0, p_D = 0
\]
условие на бесконечном расстоянии возмущения от скважине нет \[
r_D = \infty, p_D = 0
\]
В этом случае решение может быть выражено через функцию интегральной экспоненты \[
p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right)
\]
где -Ei(-x) - интегральная показательная функция.
Решение в размерных переменных можно записать как \[
p\left(r,t\right)=p_i-\frac{18.41q_sB\mu}{kh}\left(-\frac{1}{2} Ei \left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right)
\]
import scipy.special as sc"""Решение линейного стока уравнения фильтрации"""def pd_ei(td, rd):""" Решение линейного стока уравнения фильтрации rd - безразмерное расстояние td - безразмерное время """return-1/2*sc.expi(-rd**2/4/ td)
Рисунок 16: Решение линейного стока в безразмерных переменных
5.1 Построение графиков решения в размерных координатах
Для расчетов в размерных переменных необходимо определить функции перевода размерных координаты в безразмерные и обратно. В целом расчеты станут более громоздкие.
""" определим функции для перевода размерных переменных в безразмерные и обратно пригодится для построения графиков и ведения расчетов при наименовании функций придерживаемся следующих соглашений сначала идет название того, что считаем в конце указывается размерность результата, если это уместно"""def r_from_rd_m(rd, rw_m=0.1):""" перевод безразмерного расстояния в размерное rd - безразмерное расстояние rw_m - радиус скважины, м """return rd*rw_mdef rd_from_r(r_m, rw_m=0.1):""" перевод размерного расстояния в безразмерное r_m - размерное расстояние, м rw_m - радиус скважины, м """return r_m/rw_mdef t_from_td_hr(td, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):""" перевод безразмерного времени в размерное, результат в часах td - безразмерное время k_mD - проницаемость пласта, мД phi - пористость, доли единиц mu_cP - динамическая вязкость флюида, сП ct_1atm - общая сжимаемость, 1/атм rw_m - радиус скважины, м """return td * phi * mu_cP * ct_1atm * rw_m * rw_m / k_mD /0.00036def td_from_t(t_hr, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):""" перевод размерного времени в безразмерное t_hr - размерное время, час k_mD - проницаемость пласта, мД phi - пористость, доли единиц mu_cP - динамическая вязкость флюида, сП ct_1atm - общая сжимаемость, 1/атм rw_m - радиус скважины, м """return0.00036* t_hr * k_mD / (phi * mu_cP * ct_1atm * rw_m * rw_m) def p_from_pd_atma(pd, k_mD=10, h_m=10, q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):""" перевод безразмерного давления в размерное, результат в абсолютных атмосферах pd - безразмерное давление k_mD - проницаемость пласта, мД h_m - мощность пласта, м q_sm3day - дебит на поверхности, м3/сут в с.у. fvf_m3m3 - объемный коэффициент нефти, м3/м3 mu_cP - динамическая вязкость флюида, сП pi_atma - начальное давление, абсолютные атм. """return pi_atma - pd *18.41* q_sm3day * b_m3m3 * mu_cP / k_mD / h_m def pd_from_p(p_atma, k_mD=10, h_m=10, q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):""" перевод размерного давления в безразмерное p_atma - давление k_mD - проницаемость пласта, мД h_m - мощность пласта, м q_sm3day - дебит на поверхности, м3/сут в с.у. fvf_m3m3 - объемный коэффициент нефти, м3/м3 mu_cP - динамическая вязкость флюида, сП pi_atma - начальное давление, абсолютные атм. """return (pi_atma - p_atma) / (18.41*q_sm3day*b_m3m3*mu_cP) * k_mD * h_m
решение с интегральной экспонентой можно задать используя scipy.special или сгенерировать из решения sympy
pd_ei_ = sp.lambdify([t_d, r_d],eq1.rhs, modules=["scipy","numpy"]) # сравним решения заданные явно и сгенерированные из решения `sympy`print(pd_ei(1000,1))print(pd_ei_(1000,1))print(pd_ei(1000,1)==pd_ei_(1000,1))
Решение в размерных переменных можно записать как \[
p\left(r,t\right)=p_i-\frac{18.41q_sB\mu}{kh}\left(-\frac{1}{2} Ei \left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)\right)
\]
Для совершенствования навыков работы с python выполните следующие задания:
Постройте графики зависимости забойного давления от времени в полулогарифмических координатах
Сравните графики распределения давления вокруг скважины с использованием стационарного решения и решения линейного стока
Постройте график распределения давления в пласте для композитного пласта. В композитном пласте на расстоянии \(r<r_1\) проницаемость равна \(k=k_1\), а для \(r>=r_1\), \(k=k_2\). Как будет меняться воронка депрессии в таком пласте со временем?
Постройте двумерную тепловую карту распределения давления в пласте для моделей однородного пласта и композитного пласта.
Постройте график зависимости радиуса влияния скважины от времени.
Постройте решение для линейного источника на плоскости проинтегрировав точечное решение по координате. Постройте поле давления для такого источника. Считайте дебит линейного источника известным.
6 Расчет кривой восстановления давления
Один из самых простых примеров применения суперпозиции. Предполагаем, что добывающая скважина в однородном изотропном пласте запускается в момент времени t=0 и работает t_p_hr часов, после чего останавливается. После остановки скважины забойное давление растет - и мы получим кривую восстановления давления.
Пусть решение задачи запуска скважины (падения давления) будет \(P_D(t_D, r_D)\). Тогда решение для изменения давления при запуске и последующей остановки скважины можно представить в виде \[\begin{equation}
P_{bu.D}(t_D, t_{prod.D}, r_D) = P_D(t_D) - P_D(t_D-t_{prod.D}, r_D) \cdot \mathcal{H}(t_D-t_{prod.D})
\end{equation}\]
где
\(t_D\) - безразмерное время после запуска скважины,
\(t_{prod.D}\) - безразмерное время работы скважины после запуска
\(\mathcal{H}\) - ступенчатая функция Хевисайда (в некоторых книгах обозначается как \(\theta\))
\(P_D(t_D, r_D)\) - безразмерное давление - решение задачи запуска скважины (падения давления)
\(P_{bu.D}(t_D, t_{prod.D}, r_D)\) - безразмерное давление- решение задачи запуска скважины и последующей остановки скважины
Для проведения векторных расчетов в python удобно выражение с использованием функции Хевисайда
\[ \mathcal{H} = \begin{cases}0 & x < 0\\1 & x = 0\\1 & x > 0\end{cases}\]
Применение функции Хевисайда позволяет избежать в расчетных функциях применение условных операторов в явном виде для отдельных элементов входных массивов. Это потенциально ускоряет расчет.
def pd_build_up(td, td_p, rd):""" расчет давления для запуска и последующей остановки скважины td - время после запуска td_p - время безразмерное - которое скважина работала до остановки rd - расстояния от скважины """# применение функции Хевисайда здесь делает расчет корректным# для входных векторов tdreturn pd_ei(td, rd) - np.heaviside(td-td_p,1) * pd_ei(td-td_p, rd)
Show the code
t_arr = np.logspace(-10, 2, 1000)t_prod_hr =24k =10# проницаемостьq =30# дебит# переведем размерный массив времени в безразмерные величины# некоторые параметры можно было бы пропустить и оставить значения по умолчанию# но для полноты приведем их в явном видеtd_arr = td_from_t(t_arr, k_mD=k, phi=0.2, mu_cP=1, ct_1atm=1e-05, rw_m=0.1)td_prod = td_from_t(t_prod_hr, k_mD=k, phi=0.2, mu_cP=1, ct_1atm=1e-05, rw_m=0.1)print('время работы скважины {:.2f} часа, что соответсвует безразмерному времени {:.2f}'.format(t_prod_hr, td_prod))# для заданного массива безразмерных времен рассчитаем безразмерные давленияpd_arr = pd_build_up(td_arr, td_prod, rd=1)# построение графикаplt.rcParams["figure.figsize"] = (8,3)plt.plot(td_arr, pd_arr)plt.xlabel('td')plt.ylabel('pd')plt.title('Решение в безразмерных координатах')plt.show()
время работы скважины 24.00 часа, что соответсвует безразмерному времени 4320000.00
Рисунок 21: График изменения забойного давления при запуске и остановке скважины в безразмерных координатах
Рисунок 25: График изменения давления в пласте при запуске и остановке скважины
7 Решение для произвольной истории дебитов (ступенчатое изменение дебита)
Для расчета изменения давления при переменном дебите введем произвольное референсное значение дебита $ q_{ref} $ (например первое не нулевое значение дебита при запуске скважины). Используем это значение для определения безразмерного давления. \[
p_D = \frac{kh}{ 18.41 q_{ref} B \mu} \left( p_i - p \right)
\]
и безразмерного дебита
\[
q_D = \frac{q}{q_{ref}}
\]
Тогда, используя принцип суперпозиции, можем выписать выражение для изменения давления на скважине и вокруг нее для произвольного момента времени
\(i\) - индекс значения дебита в таблице изменения дебитов
\(q_{D(i)}\) - безразмерный дебит с номером \(i\), который стартует в момент времени \(t_i\). Для первого момента времени \(i\) дебит следующий перед ним считается равным нулю
\(t_{D(i)}\) - безразмерный момент времени - включения дебита с номером \(i\)
\(t_{D}\) - безразмерный момент времени для которого проводится расчет
\(p_D\left(t\right)\) - зависимость безразмерного давление от времени - решение задачи запуска скважины с постоянным единичным дебитом
$P_{mr.D} $ - безразмерное давление \(P_{mr.D}(t_D, r_D)\) учитывающее историю изменения дебитов скважины
# создадим историю изменения дебитов t_history = np.array([ 0., 2., 24. ], dtype=np.float64) q_history = np.array([10., 5., 0.], dtype=np.float64)# массивы должны быть одной и той же длиныdef pd_superposition(td, td_hist, qd_hist, rd):""" расчет безразмерного давления для последовательности безразмерных дебитов td - время расчета после запуска, безразмерное td_hist - массив времен изменения режимов работы скважин, безразмерное qd_hist - массив дебитов установленных после изменения режима работы, безразмерное """# принудительно добавим нули во входные массивы, # чтобы учесть запуск скважины qdh = np.hstack([0, qd_hist]) tdh = np.hstack([0, td_hist])# построим дебиты виртуальных скважин - # разности реальных дебитов при переключении delta_qd = np.hstack([0, np.diff(qdh)])# векторная магия - время как вектор и переключения дебитов тоже # организуем сумму по временам, внутри сумма по переключениям# делаем при помощи расчета meshgrid и поиска накопленных сумм qd_v, td_v =np.meshgrid(delta_qd, td)# используем куммулятивную сумму numpy для того что суммировать dpd = np.cumsum(qd_v * pd_ei((td_v - tdh), rd) * np.heaviside((td_v - tdh), 1), 1)# последний столбец - полная сумма, которая нужна в качестве результатаreturn dpd[:,-1]def q_superposition(t, t_hist, q_hist):""" расчет давления для запуска и последующей остановки скважины t_hr - время после запуска в часах t_hist_hr - массив времен изменения режимов работы скважин q_hist_sm3day - массив дебитов установленных после изменения режима работы """def interpolate_constant(x, xp, yp): indices = np.searchsorted(xp, x, side='right') y = np.concatenate(([0], yp))return y[indices] q=[]for ti in t: q.append(interpolate_constant(ti, t_hist, q_hist))return q
Для линейно меняющегося дебита во времени (как и для любой другой зависимости) надо решение проинтегрировать по времени (надо бы подробнее расписать - сделать это позже, например как у Щелкачева в основах нестационарной фильтрации на стр 321).
Для линейной зависимости дебита от времени \[
Q_D = dQ_D \cdot t_D
\] можно получить выражение \[
p_D(r_D,t_D, dQ_D) =-\frac{dQ_D t_D }{2} \left[ \left( 1+ \frac{r_D^2}{4 t_D} \right) Ei \left(- \dfrac{r_D^2}{4t_D} \right) + e^{-\dfrac{r_D^2}{4t_D}} \right]
\]
где \(dQ_D\) - скорость изменения дебита.
Для таблично заданных дебитов и времен можно оценить
\(p_D\left(t\right)\) - зависимость безразмерного давление от времени - решение задачи запуска скважины с постоянным единичным дебитом
$P_{mr.D} $ - безразмерное давление \(P_{mr.D}(t_D, r_D)\) учитывающее историю изменения дебитов скважины
следует обратить внимание, при суперпозиции скорость изменения дебита вычисляется как \(dQ_{D(i+1)} - dQ_{D(i)}\). при реализации расчета необходимо предусмотреть, чтобы для первого и последнего шага расчет прошел корректно. Для этого можно, например, добавить к массивам дебитов и времени дополнительный значения в начале и в конце массивов соответствующие постоянным значениям дебита.
Также надо учитывать, что в приведенном выражении массивы должны начинаться со значений \(Q_D=0\)
# Решение линейного стока уравнения фильтрацииdef pd_lin_ei(td, rd=1, dqd_dtd=1):""" Решение линейного стока уравнения фильтрации rd - безразмерное расстояние td - безразмерное время """# при расчете убедимся, что td=0 не повлияет на расчет, # даже если td массив и нулевой только один элемент td = np.array(td, dtype =float) pd = (1+ rd**2/4/td) * (-sc.expi(-rd**2/4/td)) - np.exp(-rd**2/4/td)return dqd_dtd * td * pd /2def pd_superposition_lin(td, td_hist, qd_hist):""" расчет безразмерного давления для последовательности безразмерных дебитов - td - время расчета после запуска, безразмерное - td_hist - массив времен изменения режимов работы, безразмерное - qd_hist - массив дебитов после изменения режима работы, безразмерное """# принудительно добавим нули во входные массивы, # чтобы учесть запуск скважины qdh = np.hstack([qd_hist]) tdh = np.hstack([td_hist])# построим дебиты виртуальных скважин # разности реальных дебитов при переключении delta_qd = np.hstack([np.diff(qdh),0]) delta_td = np.hstack([np.diff(tdh),1]) dq_dt = delta_qd / delta_td dq_dt = np.diff(np.hstack([0, delta_qd / delta_td]))# референсный безразмерный дебит это 1# векторная магия - время как вектор и переключения дебитов тоже # организуем сумму по временам, внутри сумма по переключениям# делаем при помощи расчета meshgrid и поиска накопленных сумм qd_v, td_v =np.meshgrid(delta_qd, td) dpd = np.cumsum(pd_lin_ei((td_v - tdh), dqd_dtd=dq_dt) * np.heaviside((td_v - tdh), 1),1 )return dpd[:,-1]def p_superposition_lin_atma(t_hr, t_hist_hr, q_hist_sm3day, k_mD=10, h_m=10, b_m3m3=1.2, mu_cP=1, pi_atma=250, phi=0.2, ct_1atm=1e-05, rw_m=0.1):""" расчет давления для запуска и последующей остановки скважины t_hr - время после запуска в часах t_hist_hr - массив времен изменения режимов работы скважин q_hist_sm3day - массив дебитов установленных после изменения режима работы k_mD=10 - проницаемость, мД, h_m=10 - мощность пласта, м, b_m3m3=1.2 - объемный коэффициент, м3/м3, mu_cP=1 - вязкость нефти, сП, pi_atma=250 - начальное давление, атм, phi=0.2 - пористость, доли единиц, ct_1atm=1e-05 - общая сжимаемость, 1/атм, rw_m=0.1 - радиус скважины """ q_ref =1. td = td_from_t(t_hr, k_mD=k_mD, phi=phi, mu_cP=mu_cP, ct_1atm=ct_1atm, rw_m=rw_m) td_hist = td_from_t(t_hist_hr, k_mD=k_mD, phi=phi, mu_cP=mu_cP, ct_1atm=ct_1atm, rw_m=rw_m)return p_from_pd_atma(pd_superposition_lin(td, td_hist, q_hist_sm3day / q_ref), k_mD=10, h_m=10, q_sm3day=q_ref, b_m3m3=1.2, mu_cP=1, pi_atma=250)
Рисунок 30: График изменения давления в пласте при работе скважины с разными дебитами
8 Радиус влияния скважины
Нестационарное решение в безразмерных переменных \[
p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right)
\] где безразмерные переменные введены как \[
r_D = \frac{r}{r_w}
\]\[
t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2}
\]\[
p_D = \frac{kh}{ 18.41 q_s B \mu} \left( p_i - p \right)
\]
Здесь использование единицы измерения СИ. - \(r_w\) - радиус скважины, м - \(r\) - расстояние от центра скважины до точки в пласте, м - \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям м3/сут - \(\phi\) - пористость, доли единиц - \(\mu\) - вязкость нефти в пласте, сП - \(B\) - объемный коэффициент нефти, м3/м3 - \(p_i\) - начальное давление в пласте, атм - \(p\) - давление на расстоянии \(r\), атм - \(c_t\) - общая сжимаемость системы в пласте, 1/атм
Для этих же безразмерных переменных, считая начальное давление равным давлению на контуре можно записать стационарное решение для движения в круговом пласте
\[
p_D = \ln r_{eD} - \ln r_D
\]
сравним это решение с логарифмической аппроксимацией (1)
которое можно преобразовать к виду \[
p_D(r_D,t_D) = - q_D \ln r_D + \frac{q_D}{2} \left[ \ln(4t_D) -\gamma \right]
\]
сравнивая со стационарным решением можно найти выражение безразмерного радиуса контура в зависимости от безразмерного времени \[
\ln r_{eD} = \frac{1}{2}(\ln(4t_D)-\gamma)
\]
\[
r_{eD} = \sqrt { 4t_D e^{-\gamma} }
\]
# оценим значение величины под корнемprint(4*np.exp(-0.57721566481))
2.2458379344731085
наконец получим \[
r_{eD} = \sqrt {2.2458 t_D}
\]
это значение называют радиусом влияния скважины. Используя это значение для определенного момента времени можно получить стационарное распределение давления в системе хорошо приближающее решение линейного стока работающего в бесконечном пласте. Можно считать это расстояние за расстояние на которое распространяется влияние скважины.
достижение радиуса влияния внешних границ будет обуславливать начало перехода от неустановившегося режима фильтрации к режиму обусловленному влиянием границ - стационарному для границы постоянного давления или псевдоустановившемуся для границы непротекания.
# Решение линейного стока уравнения фильтрацииdef pd_ei(td, rd):""" Решение линейного стока уравнения фильтрации - rd - безразмерное расстояние - td - безразмерное время """return-1/2*sc.expi(-rd**2/4/ td)def pd_ss(rd, red):""" стационарное решение в безразмерных переменных - rd - безразмерное расстояние - red - безразмерное расстояние до границы """return np.log(red/rd)
Show the code
# зададим точки расстояний для отрисовки графикаrdl = np.logspace(-1, 3 , 100)fig, (ax1, ax2) = plt.subplots(1, 2, figsize = [12,5])# построим первый график в обычных координатахfor td in [10, 100, 1000]: ax1.plot(rdl, pd_ei(td, rdl)) red = (td*2.2458)**0.5 ax1.plot(rdl, pd_ss(rdl, red))# построим второй график в полулогарифмических координатахfor td in [10, 100, 1000]: ax2.plot(rdl, pd_ei(td, rdl)) red = (td*2.2458)**0.5 ax2.plot(rdl, pd_ss(rdl, red))plt.xscale('log')plt.show()
Рисунок 31: Про радиус влияния скважины
9 Преобразования размерных величин с использованием python и sympy
Преобразования размерных величин удобно выполнять с модулем символьных вычислений python - sympy. Преобразования размерностей ключевых величин полезно знать наизусть, хотя всегда можно найти их в таблицах. Значения многих физические константы зашины в модуле scipy.constants, иногда это оказывается удобным, при этом автоматически будет учитываться достаточно большое количество знаков после запятой в константах. Рассмотрим размерности ряда величин широко применяемых в нефтяном инжиниринге.
# импортируем sympy и константы из scipyimport sympy as spimport scipy.constants as const# в модуле scipy.constants есть значения общепринятых констант -например значение piprint(const.pi)
3.141592653589793
9.1 Задание размерных величин
9.1.1 Объемный расход \(q\)
В СИ измеряется в [м\(^3\)/сек], в практических метрических единицах измеряется в [м\(^3\)/сут], в американских промысловых единицах измеряется в [bbl/day].
# выведем некоторые переводные коэффициенты для объемных расходовprint(f'Одни [сут] = {24*60*60} = {const.day} [сек]')print(f'Один [м3/сут] = {1/const.day} [м3/сек]')print(f'Один баррель в день [bbl/day] = {const.bbl} [м3/сут]')print(f'Один баррель в день [bbl/day] = {const.bbl/const.day} [м3/сек]')print(f'Один [м3/сут] = {1/const.day} [м3/сек]')
Одни [сут] = 86400 = 86400.0 [сек]
Один [м3/сут] = 1.1574074074074073e-05 [м3/сек]
Один баррель в день [bbl/day] = 0.15898729492799998 [м3/сут]
Один баррель в день [bbl/day] = 1.8401307283333331e-06 [м3/сек]
Один [м3/сут] = 1.1574074074074073e-05 [м3/сек]
9.1.2 Проницаемость \(k\)
В СИ измеряется в [м\(^2\)], в практических метрических единицах измеряется в [мД], в американских промысловых единицах измеряется в [mD].
Определение: в пористой среде с проницаемостью в один Дарси для поддержания течения жидкости с динамической вязкостью 1 сП со скоростью фильтрации 1 см/с необходимо поддерживать перепад давления жидкости приблизительно в одну атмосферу на 1 см вдоль направления течения. При использовании физической атмосферы для расчета перепада давления (физическая атмосфера = 101 325 Па) единица проницаемости равняется приблизительно 0.986923 мкм².
В отечественной литературе при определении дарси в качестве величины атмосферы было принято использовать техническую атмосферу (1 кгс/см² = 98 066,5 Па), так что для величины дарси получалось значение приблизительно 1,02 мкм², причём эпизодические случаи использования западного определения дарси специально отмечались [ru.wikipedia.org/wiki/Дарси]. Согласно ГОСТ 26450.2-85 величины 1 Дарси \(= 0.9869⋅10^{−12}\) м².
AT =98066.5# technical atmosphere in Pa, техническая атмосфера в Паprint(f'Один [psi] в [Па] = {const.psi}')print(f'Один [bar] в [Па] = {const.bar}')print(f'Один [atm] в [Па] = {const.atm}')print(f'Один [at] в [Па] = {AT}')print(f'Один [atm] в [psi] = {const.atm/const.psi}')
Один [psi] в [Па] = 6894.757293168361
Один [bar] в [Па] = 100000.0
Один [atm] в [Па] = 101325.0
Один [at] в [Па] = 98066.5
Один [atm] в [psi] = 14.69594877551345
9.1.5 Расстояние \(x\)
\(1\) [м] = \(3.28\) [ft]
9.2 Размерный коэффициент для формулы Дюпюи
Используя рассчитанные выше переводные коэффициенты для различных размерных величин рассчитаем переводной коэффициент в формуле Дюпюи
# зададим переменные sympyQ, k, h, mu, B, pres, pwf, re, rw, S, pi = sp.symbols('Q k h mu B p_res p_wf r_e r_w S pi', real=True, positive=True)# определим уравнениеeq = sp.Eq(Q, 2* pi * k * h / (mu * B) * (pres - pwf) / (sp.ln(re/rw) + S))display(eq)
\(\displaystyle Q = \frac{2 h k \pi \left(p_{res} - p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
eq = eq.subs(Q, 1/const.day * Q) # дебит, [м3/сут] в [м3/сек]f_k =1e5/const.atm *1e-15eq = eq.subs(k, f_k * k) # проницаемость, [мД] в [м2]eq = eq.subs(mu, 1e-3* mu) # вязкость, [сП] в [Па сек]eq = eq.subs(pres, const.atm * pres) # давление [атм] в [Па]eq = eq.subs(pwf, const.atm * pwf) # давление [атм] в [Па]eq = eq.subs(pi, const.pi)display(eq)
\(\displaystyle Q = \frac{0.0542867210540316 h k \left(p_{res} - p_{wf}\right)}{B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Выделим полученную константу в явном виде и найдем обратную величину - это и будет необходимый нам переводной коэффициент.
f =1/eq1.args[0]print(f)
18.4207110060064
По умолчанию sympy автоматически организует порядок элементов в своих выражениях. Этот порядок может отличаться от привычного - хотя и суть формул при этом не меняется. Применяя некоторые хитрости можно заставить sympy вывести выражения в приемлимом виде.
a = sp.symbols('a')eq2 = eq1.subs(eq1.args[0],1/a)with sp.evaluate(False): display(eq2.subs(a, f))
\(\displaystyle \frac{h k \left(p_{res} - p_{wf}\right)}{18.4207110060064 B \mu \left(S + \log{\left(\frac{r_{e}}{r_{w}} \right)}\right)}\)
Но иногда результат проще переписать руками в нужном виде. В итоге уравнение Дюпюи в практических метрических единицах измерения примет вид.